#include "energybands.h"


EnergyBands::EnergyBands(void) :
	// A pointer to the used system. The system contains variables 
	// needed to construct elements (coupling strengths etc.).
	system(&System::Std),
	// A double containing the current bias (for bi-layered 
	// systems).
	bias(0.),
	// A enum type containing the current k-path to use. It can be
	// either brillouinpath or userdefined.
	k_val(brillouinpath),
	// A double matrix containing the current k-points.
	k(mat()),
	// A enum type containing the what eigen values to generate.
	// This can be all, count (a spec. number of eig), or interval
	// (a spec. value interval)
	eig_val(all),
	// A size_t containing the amount of eigvalues to calculate. If
	// this count is to big it will get truncated at calc. time.
	N_eig(20),
	// Two doubles containing the interval of eigvalues to calculate.
	eig_min(-1.), eig_max(1.)
{ }

#include <mkl.h>

EnergyBands &EnergyBands::operator()(int id)
{
	
	// First we get the unit from the id.
	Unit* unitref = system->getUnit(id);
	if(!unitref) // If the unit is invalid then we do nothing.
		return *this;

	// First we get the k path.
	if(k_val == brillouinpath)
		k = getBrillouin(id);
	else if(k.empty())
		return *this;

	// We prepare the hamilton.
	Hamiltonian hamilton;
	hamilton.system = system;
	hamilton.bias = bias;

	char range;
	
	// Now we allocated for the result.
	if(eig_val == all || eig_val == interval)
	{
		range = (eig_val == all) ? 'A' : 'V';
		bands.set_size(k.n_cols,unitref->getAtoms().size());
	}
	else 
	{
		range = 'I';
		bands.set_size(k.n_cols,N_eig);
	}


	
	for(size_t nk = 0; nk < k.n_cols; nk++)
	{
		hamilton.k = k.col(nk);
		cx_mat H = cx_mat(hamilton(id).H);

		vec energies(H.n_cols);

		long long m; // The number of eig_vals found!
		long long* isuppz = new long long[2*bands.n_cols];

		long result = LAPACKE_zheevr(
			LAPACK_COL_MAJOR, 'N', range, 'U', 
			H.n_cols, (MKL_Complex16*)H.memptr(), H.n_rows, 
			eig_min, eig_min,
			(H.n_cols+1)/2-(N_eig-1)/2, 
			(H.n_cols+1)/2+(N_eig)/2,
			0., &m, energies.memptr(), 0, H.n_cols, isuppz );
		
		if(result == 0)
		{
			for(size_t en = 0; en < bands.n_cols; en++)
			{
				if(m > en)
					bands(nk,en) = std::pair<vec2,double>(k.col(nk).rows(0,1),energies(en));
				else
					bands(nk,en) = std::pair<vec2,double>(k.col(nk).rows(0,1),0.);
			}
		}
	}

	return *this;
}

mat EnergyBands::getBrillouin(int id) const
{
	mat kpath;
		
	// First we get the unit from the id.
	Unit* unitref = system->getUnit(id);

	const size_t size = unitref->getBrillouin().n_rows;

	if(unitref->getPeriodicity().n_cols == 1)
	{
		kpath.set_size(size,2);
		kpath.col(0) = unitref->getBrillouin().col(0);
		kpath.col(1) = unitref->getBrillouin().col(1);
	}
	else if(unitref->getPeriodicity().n_cols == 2)
	{
		kpath.set_size(size,4);
		kpath.col(0) = unitref->getBrillouin().col(0);										// K
		kpath.col(1) = (unitref->getBrillouin().col(0)+unitref->getBrillouin().col(1))/2;	// M
		kpath.col(2).fill(0.);																// G
		kpath.col(3) = unitref->getBrillouin().col(0);										// K
	}

	return interpolatePath(kpath,100);
}

mat EnergyBands::interpolatePath(const mat &path, const long &points)
{
	long pathsize = path.n_cols;
	mat paths(path.n_rows,pathsize-1);
	vec lengths(pathsize-1);
	uvec pointcounts(pathsize-1);
	
	for(long p = 0; p < pathsize - 1; p++)
		paths.col(p) = path.col(p+1)-path.col(p);
	for(long p = 0; p < pathsize - 1; p++)
		lengths(p) = arma::norm(paths.col(p),2);

	double fulllength = arma::sum(lengths);
	
	for(long p = 0; p < pathsize - 1; p++)
		pointcounts(p) =  ((pathsize > 2) ? points : points-1)*lengths(p)/fulllength;

	mat result(path.n_rows,points);
	size_t currentpoint = 0;

	for(size_t p = 0; p < pathsize-1 && currentpoint < points; p++)
	{
		for(size_t pp = 0; pp < pointcounts[p] && currentpoint < points; pp++)
		{
			double f = double(pp)/double(pointcounts(p));
			result.col(currentpoint) = (1.-f)*path.col(p) + f*path.col(p+1);
			currentpoint++;
		}
	}

	if(currentpoint < points)
		result.col(currentpoint) = path.col(pathsize-1);
	
	return result;
}